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Abstract 

The Principal results of a recent theory of fuel optimal space trajectories for linear differential equations are 
presented. Both impulsive and bounded-thrust problems are treated. A new form of the Lawden Primer vector 
is found that is identical for both problems. For this reason, starting iteratives from the solution of the impulsive 
problem are highly effective in the solution of the two-point boundary-value problem associated with bounded 
thrust. These results were applied to the problem of fuel optimal manuvers of a spacecraft near a satellite in circular 
orbit using the Clohessy-Wiltshire equations. For this case two-point boundary-value problems were solved using 
a microcomputer, and optimal trajectory shapes displayed. The results of this theory can also be applied if the 
satellite is in an arbitrary Keplerian orbit through the use of the Tschauner-Hempel equations. A new form of the 
solution of these equations has been found that is identical for elliptical, parabolic, and hyperbolic orbits except in 
the way that a certain integral is evaluated. For elliptical orbits this integral is evaluated through the use of the 
eccentric anomaly. An analogous evaluation is performed for hyperbolic orbits. 
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1 Introduction 


Many problems of optimal trajectories, maneuvers, and rendezvous of spacecraft have been investigated using linearized equations 
of motion. Linear equations describing the relative motion of a spacecraft near a satellite in circular orbit [1-4] have been very 
useful to aerospace researchers. Some of the published applications are cited here [5-8]. The linear equations describing relative 
motion near a satellite in an elliptical orbit [9,10] are not quite so well known, but have also been useful in applications [11- 
16]. These equations also generalize to a description of relative motion near a satellite in an arbitrary Keplerian orbit [17,18]. 
Another approach to similar problems involves the linearization with respect to the orbital parameters associated with a satellite 
or spacecraft [19]. Other sets of linearized equations describing the motion of an object near one of the five Lagrange points in 
the restricted three body problem are well known in celestial mechanics [20]. 

In this paper we present necessary and sufficient conditions for fuel-optimal trajectories of a spacecraft whenever the 
equations of motion are linear. This theory encompasses all of the preceeding examples in which the equations of motion of a 

spacecraft are usefully linearized about an equilibrium condition [1-20], 

The work divides naturally into two distinct areas, impulsive problems, and bounded thrust problems, and is based largely 
on recent investigations [21,22] in these two areas. Early work on these problems was done by Neustadt [23,24] who formulated 
the linear impulsive spacecraft trajectory problem mathematically as an unbounded thrust problem in a class of more general 
nonlinear programming problems. He obtained an existence theorem and necessary conditions for solution of this problem. He 
presented also a precise sense in which the impulsive problem solution is a limit of the bounded thrust problem solution as the 
bound becomes arbitrarily large, and showed that the necessary conditions for the impulsive problem are obtained from the 
necessary conditions for the bounded thrust problem by passing to the limit. 

Our approach differs significantly from Neustadt’s, especially in the impulsive problem. We formulate this problem through 
the use of a finite number of velocity increments as independent variables instead of thrust functions having unbounded range. 
This results in a simpler problem that can be solved without the use of mathematical control theory or advanced mathematics. 
The bounded thrust problem is then solved by the well known Principle of Pontryagin. It is found that the two problem solutions 
are closely related, and insight into either problem is sometimes obtained from the other. Probably the most useful relationship 
is the fact that both problems are found to contain the identical new form of the primer vector function originally defined by 
Lawden [25]. 

It is shown that a class of two-point boundary-value problems associated with bounded thrust can be readily solved from 
starting iteratives that are obtained from a primer vector function associated with the impulsive problem. 

Section two demonstrates necessary and sufficient conditions for solution of the impulsive minimum fuel problem with 
linear equations of motion. At the end of the section this material is applied to the case of a spacecraft near a circular orbit 
using the Clohessy- Wiltshire [1-4] equations of motion. Simulations of impulsive spacecraft trajectories are presented. In section 
three, necessary and sufficient conditions are also revealed for solution of the related bounded thrust minimum fuel problem. 
Simulations of optimal bounded thrust spacecraft trajectories near circular orbit are also presented for this problem using the 
same equations of motion. The final section indicates how this material can be applied for a spacecraft near a non circular 
Keplerian orbit. The Tschauner-Hempel equations of motion [9,10] are applicable to this problem. The work of Tschauner and 
Hempel and Weiss [15] has demonstrated that solutions of these equations can be used to define a fundamental matrix solution 
that applies to either circular or elliptical orbits. This work was found useful in constructing two-impulse solutions to rendezvous 
problems involving objects in elliptical orbits of high eccentricity [16] and can also be used to apply the theory presented here 
for rendezvous of a spacecraft with objects in circular or elliptical orbits. We complement this work by presenting a new form 
[18] of the solution of the Tschauner-Hempel equations that is valid for elliptical, parabolic, or hyperbolic orbits, and that avoids 
the removable singularities found in some earlier papers. This form of solution is then used to construct a fundamental matrix 
solution. Based on this fundamental matrix solution, the theory of sections two and three can be applied to the problem of fuel 
optimal rendezvous of a spacecraft with an object near a satellite in a general Keplerian orbit. 

2 The Impulsive Minimization Problem 

2.1 The Impulsive Problem Formulation 

We let m and n be positive integers, and xo, i>o> , v/, Av; (t = are all elements of the m-dimensional Euclidean 

space & m . The real numbers 9 0 and 9 f define a bounded interval 0 = [0 o ,6j]. We denote the generalized instantaneous position 
and velocity of a spacecraft respectively by x(0) and v(0) in where 9 £ 0 represents a convenient independent variable 
such as the normalized time in flight, applications the true anomaly of a satellite in Keplerian orbit. [10,17] A prime is used 
to indicate differentiation with respect to 9 . The norm or magnitude of any vector will be represented by the symbol | • | • 
The symbol x will be used with reference to rows and columns of a matrix or will denote the Cartesian Product of sets. The 
superscript T refers to the transpose of a matrix or vector. For each 9 G 0 the vector y(9) £ refers to the transpose^ of the 
pair (x(0) T , v(0) T ). Similarly the vectors yo and yf in SR 2 ”* denote respectively the transposes of (x 0 ,t>o ) an< ^ ( x /» v / )* We 
assume that a velocity increment A Vi is added to the velocity v(0;) of the spacecraft at a point 0, (* = 1> •••> n )* This assum ption 
represents an idealization of the effect of the application of n short duration thrusts to the spacecraft, resulting in instantaneous 
changes in velocity without change in position. This type of idealization is in common practice by mission planners and can 
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be useful in providing approximate data for iterative schemes to solve problems having more accurate models. The n points of 
application of the velocity increments can be specified a priori by the mission planner. Another possibility is that the locations 

of these velocity increments be determined through the optimization process. Our approach is flexible enough to include either 
of these possibilities or a combination of the two. 

^Specifically we let p be a positive integer that represents the number of elements in 0 that are specified a priori, and 

n p - Wi G © p is this ordered set of specified points in © where velocity increments are to be applied. For technical 

reasons we shall require that both end points 9 0 and 9, are in N p so that p > 2. We also specify an integer r > 0, and we shall 
seek an ordered set N r - (9 kl , ..., 0 kr ) g 0 r of r elements of © that are to be determined through the optimization process, and 
that represent optimum locations for the application of velocity increments. If r = 0 then N T is empty. We now set n = p + r 
order the components of (N p ,N r ) € © n to form the ordered set (0,,. ,.,0„) g ©" where 9, < 9 i+1 (i = i). We observe 

that it is possible to have 0, = 0 } for » jt j. We associate a velocity increment At), with each 0, (t = 1, ...,n) even though the 
optimization process may show some of these velocity increments to be zero. Although we associate velocity increments with 
the end points 0 O and 9,, there are problems where one or both of these velocity increments become zero. We now define the 
impulsive optimal spacecraft trajectory problem for linear equations of motion as follows: 

2.1.1 Statement of the Linear Impulsive Minimization Problem 

Having specified N p and the integer r, find N r = (8 kl , ..., 9 kr ) g © r and the velocity increments At), g s )? m (i = i n ) 
which are added to the velocity vectors t )(0,) € ST* (« = 1, ...,n) in order to minimize the total characteristic velocity 

n 

J = ^|At),| (!) 

1 = 1 

of a spacecraft, whose motion is defined by the linear differential equation 

y'(0)=A(0)y(0) (2) 

for each 0 g 0 except at 9 { (» = 1, .... n), where A is a 2m x 2m real matrix-valued function continuous on ©, and whose initial 
and terminal conditions are given by 

y(0o) = y 0 , y(9j) = y f . ( 3 ) 


2.1.2 Fundamental Matrix Solution 

It is known from the theory of linear differential equations that a 2m x 2m fundamental matrix solution which we shall denote 
by $(0) is associated with the linear differential equation (2), its elements are continuously differentiable, its inverse $(0) -1 
exists and its elements are also continuously difFerentiable, for each 9 g 0. As we shall show, the 2m x m matrix which consists 
of the last m columns of 4>(0) (i. e „ the right hand half of $(0)-‘) plays an important role in the determination of the optimal 

impulses and will be denoted by R(9). Because of the continuity of A{9) the 2m x m matrix R(0) of real-valued functions is 
continuously differentiable on ©. Moreover if the elements of A{9) are analytic on ©, then the elements of R(0) are known to 
be analytic on © also. Our principal result, as presented in the next section, is stated in terms of the matrix R(9), and the first 
step in solving actual problems is to determine the matrix R($). 


2.2 Necessary and Sufficient Conditions 

As we shall see, the effects of the initial and terminal positions and velocities (3) in defining an optimal solution are determined 
exclusively through the vector b g 3i 2m which is defined as follows: 

6 = $(0 / )- 1 y / -$(0o)- 1 yo. (4) 

A family of initial and final positions and velocities having the same value of the vector b is associated with a family of problems 
having identical optimal impulsive solutions. 

In the following we defined, =| At), | (« = 1 n) and for each i such that a, / Owe define u, = At), /a,. These are simply 

the magnitudes and the unit direction vectors of the nonzero velocity increments. We now state our principal result as follows: 
Theorem 2.1 Suppose that b is nonzero, A{9) exists and is continuous on 0, and N p contains both of the end points 0 O 

Tt °l The " f° reach r = °> 1 > 2 > - andn = p+r, the element (0*, , ..., 9 kr , At>j , ..., At)„) solves the minimizatidn problem 
defined by (1 — 3) subject to the imposition of the n velocity increments, over the set © r x 9i mn if and only if there exists a 
nonzero vector A € 9i 2m such that 

cm =0 or u , = —R($i ) T A (» = 1, ..., n) 

<»i = 0or X T R(9,)R(9 i ) T X = 1 (, = l, ...,„) 

<**, = 0, or 9 kl = 0 O or 0 ki = 9, or X T R'(0 ki )R(0 k .) T X = 0 (» = 1, ..., r) 
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( 5 ) 

( 6 ) 
( 7 ) 



£[iW)tf(0i) T *]«. = ~ b 

(8) 

«* > o (» = i, ..,«) 

(9) 

= -b T X>0 

(10) 

i= 1 

b T A is a minimum on the se<{A E 3R 2m | (6 “ 10) 15 valid}. 

(11) 


A proof of this theorem is presented in [21]. The requirement that b is nonzero in this theorem prevents the problem from 
becoming trivial. If b = 0 then the minimization problem can be solved by setting At), _ 0 (t — In applying this 

theorem to actual problems one specifies the number r of points of application of velocity increments to be determined by the 
optimization process. Of course, the choice of the number r affects the solution of the minimization problem defined by (1-3) 
and the resulting minimum value of J. Although the proof introduces some additional mathematical framework and will not Tie 
considered here, it can be shown that if the elements of A(9) are analytic functions of 9 then there is an integer r* such that 
the solution of the minimization problem for r > r* does not lead to a lower value of J than does the solution for r*. One 
consequence of this fact is that N p has no effect on the minimum value of J if the integer r is sufficiently large. For this reason, if 
r is sufficiently large, then there is no loss in generality in making the assumption that N p contains the elements 9 o and 9 } . lhe 
determination of the smallest such integer r* for a wide class of problems defined by (1-3) is a problem whose solution wou 
reveal to the mission planner the maximum number of thrust impulses needed for a mission. Neustadt’s work [23] would suggest 
that r* < 2m. We can sometimes keep the value of r general. If in the determination of the locations for velocity increments, 
we obtain 9i = 9, for i / j (i.e. repeated roots) then the value of r can be lowered without altering the minimum value of J. 


2.2.1 The Primer Vector Function 

Some of the results of this theorem can be interpreted in terms of the well known primer vector theory, introduced by Lawden [25] 
and developed for impulsive spacecraft trajectory applications by Lion and Handelsman [26,27]. We shall refer to the function 
q . @ ► g£ m defined by q(0) = R(0) T X as the primer vector function. In terms of this function, (6) and (7) respectively become 

oti = 0 or | q(9,) |= 1 (i = 1 n) ( 12 ) 

o fc , = 0, or 9 ki =9 0 , or 6 k , - 6 } , or \q (9 ki ) |'= 0 (« = 1, ..., r). (13) 

Geometrically, this says that nonzero optimal impulses must occur where the magnitude of the primer is unity and, if they are 
at interior points of 6, where the magnitude of the primer is tangent to a horizontal line one unit above the horizontal axis. 
This and other related geometric conditions were observed by Lion and Handelsman, even though the class of problems under 
their investigation was very different from the above. 

2.2.2 A Priori Specification of Velocity Increments 

As previously indicated, this theorem is presented in such a way that the mission planner can specify some of the locations of the 
velocity increments a priori. It can be made even more flexible by allowing the mission planner also to specify completely some 
of the velocity increments a priori, as well. This is accomplished as follows. We assume that an additional s velocity increments 
are completely specified by the mission planner where s is a positive integer, and N p is extended to include the locations of these 
velocity increments. The total number of velocity increments is then n + s, and the s completely specified velocity increments will 

be denoted At) n+1 Av n+S with respective locations 0 n+1 ,...,0 n+s € N p . The Theorem can be restated with this adjustment 

and only minor changes in the proof are necessary. The principal change occurs in Lemma 4.1 of the next section where the 
vector b is adjusted as follows: 

b = 1 — <f>(0o) 1 yp - 

;= i 

This definition of b replaces (4) only in the following adjusted result: 

Theorem 2.2 Suppose that b , as defined above , is nonzero , A(0) exists and is continuous on 0, and that s velocity 
increments Au„+i ,...,Av„ +a E 9R m are specified a priori and N p contains their respective locations 0»+i, ..., 0 n + s and the end 
points Bo and 0 f of 0. Then for each r = 0,1,2,... and n = p - s + r, the element (0 kl , • * • , 0 kr , At>i, Ai> n ) solves the 
minimization problem defined by (1-3), subject to the imposition of the n + s velocity increments, over the set © r X « */ and 

only if there exists a nonzero vector A E such that (5-11) are satisfied where b is defined above . 
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2.3 Applications 

This approach to the solution of the minimization problem defined by (1-3) requires solving the system (5-11). The primary 
computational problem is in the solution of (6-8) which are quadratic in A and linear in (i = 1, ...,«). In most application 
to spacecraft trajectories, the position and velocity vectors of the spacecraft are in three dimensional Euclidean space so that 
m _ 3, or if the motion is restricted to an orbital plane, m = 2. If the equations of motion of the spacecraft are based on 
linearization about a satellite in circular or Keplerian orbit, we conjecture that it is not necessary to employ more than four 
impulses for one orbital period or less (0, <0 O + 2jt), so we recommend setting n = 4 for these problems. This is in agreement 
with a similar conjecture which was made for bounded thrust problems. [28,29] Neustadt’s results [23] for the general linear 
unbounded thrust problem show that the maximum number of impulses is 6 for m = 3 and 4 for planar problems. 

2.3.1 Types of Problems 

We list the following types of problems: 

1. Planar problem with four impulses with specified locations. Here (6) presents four quadratic equations in four unknowns 
to determine A. For some cases these can be solved by hand or by a computer programmed for symbol manipulation 
Half of these solutions are eliminated by the test -b T A > 0 of (10). The remaining solutions are used in (8) to solve for 
ari,...,ar n by Gauss-Jordan elimination. All that do not satisfy (9) are thrown out. Usually only one solution remains 
after this test. If more than one solution remains, the value or values of A that satisfy (11) are retained, and from this the 
velocity increment directions are determined through (5). An outstanding feature of this problem is that a change in the 
initial or terminal conditions does not require a solution of (6) again. The solutions of (6) are a set of vectors like spokes 
on a wheel, any one of which may point out the correct solution. A change in the boundary conditions defines another 
vector from the set which satisfies condition (6) and defines the optimal impulses only through the solution of sets of linear 
equations (8) to determine ori, .... or„ and then applies the tests (9-11). An example of this type of problem is presented 
in the next section. 

2. Four impulse planar problem with intermediate locations determined optimally. For this problem (6) and (7) determine 
a set of six equations in the six unknowns A, 0 2 and 0 3 . After obtaining multiple solutions of these, the tests (9) and (10) 
and the linear equations (8) can be used as in the preceding. For practical problems it may be preferable computationally 
only to approximate the solution of 0 2 and 0 3 . We guess the locations 0 2 and 0 3 and solve the preceding problem 1). We 
repeat this process several times with new values of 0 2 and 0 3 . We then pick the set which most closely approximates (7) 

. In this manner (7) can be approximated to any desired accuracy. An alternative is to use gradient methods on 0 2 and 
03 to minimize J. 

3. Four impulse three dimensional problem with specified locations of impulses. In this case (6) and (8) define a set of 10 
equations m the 10 unknowns A,ai,« 2 ,» 3 ,<* 4 . Each solution of this system is subjected to the tests (9) and (10) as in 
the previous problems. In case of multiple solutions, the test (11) is applied. 

4. A two impulse solution with impulses at the ends. Rather than using the equations (5-11) for this trivial problem, it is 
easier to use Lemma 4.1 from a recent work [21] and solve for the two impulses through the linear equations 

R(0 o )Av o + R(9f)Avf = b. ( 14 ) 

This method is computationally simpler than the two impulse method of Weiss [15] used for high eccentricity elliptic orbits. 


2.3.2 Spacecraft Maneuvers near Circular Orbit 


As an application, we consider the problem of determining optimal maneuvers of a spacecraft near a satellite in circular orbit. 
The linearized equations of motion with reference to a coordinate frame fixed in the satellite were determined independently 
by Wheelon, Clohessy and Wiltshire, Geyling, and Spradlin [1-4]. With these equations, optimal impulsive solutions were 
investigated by Prussing [6] in 1969 using primer vector theory. Our approach is to find and invert a fundamental matrix 
solution, and apply Theorem 2.1. The equations of motion when put in the state vector form yi(t) = A(t)y(t), are determined 
by the 6 x 6 matrix 


A(t) = 


0 0 
0 0 
0 0 
0 0 
0 3 w 2 
0 0 


0 1 0 0 ' 

0 0 10 

0 0 0 1 

0 0 2o» 0 ’ 

0 -2u> 0 0 

-u> 2 0 0 0 


(15) 
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Table 1: The sixteen solutions of Eq.(6) for a planar circular orbit 




1 

2 

3 

4 

1 

Ai 

0.065987656 

0.106065818 

-0.065987656 

-0.106065818 


a 2 

-0.621919006 

-0.999647 

0.621919006 

0.999647 


A3 

-0.651106246 

0.2387088 

0.651106246 

-0.2387088 


a 4 

0 

0 

0 

0 

2 

Ai 

-0.098201378 

-0.098201378 

-0.098201378 

0.098201378 


a 2 

-1.798123431 

0.052928935 

1.798123431 

-0.052928935 


A 3 

-0.182280857 

0.182280857 

0.182280857 

-0.182280857 


a 4 

0.436298624 

0.436298624 

-0.436298624 

-0.436298624 

3 

wm 




0 


m 

1 

-3 

5 

-1 

3 

5 




0 

0 



A 4 


4 

5 


^ 

11 

mm 

mmmmm 



mi 

II 

in 

I;-- Av/rVvj 

0 



1 

H 


& 

5 

_ v^E 
5 

5 

■ 

IH 


__s 

& 

& 

5 


where the real number ui is the orbital angular speed of the satellite, and the flight time is denoted by t. It is not difficult to 
find and invert a fundamental matrix $(f) associated with A(t), and obtain 


m - 1 


1 -6 0(1) 0 30(t) 2 0 

0-20 1 0 0 
0 3 sin 0(f) 0 — 2sin0(t) —cos 0(f) 0 

0 —3 cos 0(f) 0 2 cos 0(f) —sin 0(f) 0 

0 0 sin 0(f) 0 0 cos 0(f) 

0 0 cos 0(f) 0 0 -sin 0(f) 


(16) 


where 0(f) = u>t. From (16) R(t) is evident. We drop the argument f and consider 0 as the independent variable. 


We shall consider only the planar case. For this problem we find the 4x2 matrix 


m 


30 2 \ 

\ ° . 

—2 sind —cosd I 

\ 2 cosO —sinO ) 


(17) 


For simplicity we shall pick 0 o = 0,0/ = 2 tt, and we shall select four impulses at 0 ,tt/ 2 , 37 t/ 2 , and 2 x. Based on earlier studies 
for optimal trajectories with bounded thrust [ 8 ] we suspected that these locations are not far from the optimal locations. After 
computation of the A values based on these locations, this suspicion was confirmed by substituting these values of 0i into (7). 
The values of the matrix i£(0<) (« = 1 , 2 , 3 , 4 ) are easily obtained so that ( 6 ) is a manageable set of four quadratic equations 
in the four unknown components of A; Ai, A 2 , A 3 , A 4 . Solution of these four quadratic equations reveals the following 2 16 

solutions which we denote by A ij (ij = 1,2, 3 , 4 ). These solutions are presented in Table 1. We consider three distinct sets of 
initial and terminal conditions for computation of the optimal impulses using the method outlined in 1 ) of Section 2.3.1. We 
find that each of the three cases requires a different value of A from Table 1, and a different number of impulses for optimality. 
(We neglect very small impulses, and assume that they are zero because of their negligible effect in comparison with much larger 
impulses in the same case.) 

Case 1. A spacecraft is one unit behind a satellite in the same circular orbit. Its initial relative velocity is zero, and 
its object is to rendezvous with the satellite having the same velocity at 0/ = 27r. The tests (8-11) determine A12 as optimal. 
For this value of A, the solution of (8) establishes that <*2 — <*3 == 0, a two impulse solution. The optimal impulses are both 
horizontal and opposite: 

M«)= [ P | i Av(2jt) = | p ] • 

A computer simulation of the optimal trajectory is presented in Figure 1 . 

Case 2 . A spacecraft is one unit above a satellite in circular orbit with the same initial velocity as the satellite. Its object 
is to rendezvous with the satellite at 0/ — 2i r having the same final velocity. For this problem the tests ( 8 - 11 ) determined A 21 
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0.0 


Fig. 2. Optimal trajectory with three Impulses. 

X(0)T- (0.1). t*0) T « (0.0). x(2it)T . <0.0).t*2ji)T - (0.0) 

X, 



0.0 

Rg. 3. Optimal trajectory with four impluses. 
x(0) T -(1.0).xr<0)T «(0.0).x(2jt)t- (0.0). «C2n) T - (0.0.427) 
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as optimal. For this value A, eq.(8) defined a 3 = 0, a three impulse solution. The optimal impulses are 

T 1.6294 1 A T .39010 1 A _ f -.06332 1 

M°) - [ _ 66667 j - Mj) - [ .09640 J > [ -.02591 J ‘ 

A computer simulation of the optimal trajectory using these impulses is found in Figure 2. 

Case 3. A spacecraft is in circular orbit one unit behind the satellite with the same velocity as in Case 1, but in this case 
its object is to reach the satellite with a positive vertical velocity. This causes quite a different trajectory from Case 1. The 
optimal A changes to A22 and a four impulse solution results. The optimal impulses are 

, „ r -.02729 1 . ,T. r .08965 1 

M°) ~ [ .03436 J ’ A ^2^ _ [ -01194 J ’ 

3*\ r -08965 1 . ff. v r .02729 1 

My) “ [ .01194 j >M 2t ) “ [ .03436 J ' 

The optimal trajectory is presented in Figure 3. . 

Each of the figures can be compared with similar figures in [8] which are based on bounded thrust and continuous velocity. 

The similarities in the shapes are apparent. 

Extension of this work to other Keplerian orbits requires the replacement of the Clohessy-Wiltshire equations with the 
Tschauner-Hempel equations [9,10] for elliptical orbits or the equations in [18] for more general orbits. For the Tschauner-Hempel 
equations the fundamental matrix solution has been found and inverted by Weiss [15], so that the matrix R(9) can be defined. 

3 Bounded Thrust Space Trajectories 

3.1 The Bounded Thrust Minimization Problem 

We let m be a positive integer and x 0 ,v 0 ,xf,vj are elements of the Euclidean Space !R m . The real numbers t 0 and t; define a 
closed and bounded interval T = [to,*/] where t Q < t f . A dot above a variable will denote differentiation with respect to t 6 T. 
The Euclidean norm or magnitude of a vector will be referred to by the symbol | | . The superscript T denotes the transpose of 
a vector or a matrix. The elements y 0 and y f in 3? 2m are defined by yl = (lo.vj) and yj = (xf , vj). We let U denote the 
class of admissible control Junctions, which is the set of all Lebesgue measurable functions that map T into the closed unit ball 


U in 9£ m a.e. on T. 

We shall consider the linear nonhomogeneous differential equation 


y(t) = A(t)y(t) + Pw(t) 

(18) 

which is defined a.e. on T subject to the initial condition 


y(t a ) = y 0 . 

(19) 

and the terminal condition 

(20) 

y(</) = »/> 


where A is a real 2m x 2m matrix valued function on T, u»(t) T = (0 T , u(t) T ) where 0 is the zero element in 3R m ,u G U, and /? is a 
positive real number. We shall assume throughout that the elements of A are real valued analytic functions on T. A solution y 
of (18) at t G T is an element of 5R 2m , and we shall write y(t) r = (x(t) T , v(t) T ) where x(t) G and v(t) G for each t G T. 
We define the cost of any element u E U to be 

J[«]= f 1 | «(<) | dt ( 21 ) 

JtO 

which represents the Lebesgue integral of the norm of u over T. Here we are assuming that a spacecraft has constant exhaust 
velocity, and that its total fuel consumption over the interval T is proportional to J[u]. The linear bounded thrust minimization 
problem can therefore be stated as the problem of minimizing J[ti] over U subject to the conditions (18 - 20). 

Many spacecraft maneuvers and rendezvous missions can be formulated in terms of this problem. In these cases the equa- 
tions of motion of a spacecraft are approximated by the linear equations (18) where t represents a flight time or related variable 
such as the true anomaly of a satellite in Keplerian orbit, and /?u(t) usually represents an applied acceleration on the spacecraft 
caused by the engine thrust. The problem is a bounded control problem because of the restriction that u £ U (i.e. | u(t) |< 1 
a.e. on T). For simplicity, we assume the spacecraft mass is constant over T so that /? is a positive constant representing the 
maximum magnitude of the applied acceleration. We can, however, include the effect of the variation of mass as in previous 
work [29]. To avoid an additional equation that is nonlinear, we shall not do this. If the matrix A(t) is partitioned into four 
mxm matrices so that the top left is zero and the top right is the m X m identity, then i(t) = v(t) y and if t represents time in 
flight then the elements x(<) and v(t) can be interpreted respectively as relative positions and velocities of the spacecraft. 
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3.2 Necessary and Sufficient Conditions 

Before we present our principle result, we shall place a restriction on the matrix A(t) that will simplify the solution of the 

minimization problem by eliminating singular abnormal solutions. The matrix A(t) can be partitioned into four mxm matrices 
as follows: 

= ( A«(t) ^2(t) ) € T) - ( 22 ) 

We shall say that the matrix valued function A is primer-compatible if the only solution of the two equations 

A 12 (t) T p(t) = 0, p(t) = -A n (t) T p(t) (23) 

on T is the solution p(t) = 0 (t g T). It is readily seen that A is primer-compatible if Ai 2 (t) is nonsingular for each t £ T For 
this reason we see that in the usual problems in which i(t) = »(t) for each t g T the matrix A is primer- compatible because 
Ai 2 (t) — I the mxm identity matrix. By differentiating we see also that A is primer-compatible if Ai 2 (t) — vl 1 j (t)Ai 2 (t) is 
nonsingular for each t g T. Other conditions for primer-compatibility can be obtained by differentiating further. 

We now present our principal result. The proof can be found in [22]. 

Theorem 3.1 Suppose that the elements of the 2m x 2m matrix A(t) are analytic on T, and that A is primer- compatible. 
Then u is a minimum of J[u] = | «(<) | dt over U subject to the conditions (18 - SO) if and only if there exists a real number 

l o >0 and A g S? 2m such that q(t) = R(t ) T A for each t g T and 


y(t) = m 


|$(*o) -I- 0 



(fgr) 


(24) 


where $(t) is any fundamental matrix associated with A{t) for each t g T, its inverse exists and is partitioned so that $(<)“' = 
R(t)) where L(t) and R(t) are 2m x m matrices for each t g T, and either 

i) A is the zero element in » 2m ,u(<) = 0 o.e. on T, and $(t/)$(t 0 ) _1 y 0 = y/ , or else 

ii) A is nonzero, the equation | q(t) |= 0 has at most finitely many solutions on T, y(tj) = y } , and 


«(<) = - 


*(<) 

1 9 ( 0 1 


m 


a.e. on T where either the equation | q(t) |= l o /0 has at most finitely many solutions on T and 


(25) 


/(<) = 


{ 


0, I 9(0 |< to/P 

1, I ?(<)!> t o /0 


(26) 


a.e. on T, or else it is satisfied identically on T and 0 < f(t) < 1 a.e. on T. Moreover A 7 '($(t) _1 y(t) - $(f 0 ) _1 j/ 0 ) = -0 f‘ | 

9(r) | f(r)dr <0 (t g T). In particular, A T ($(t / )- 1 y / - $(t 0 ) -1 y 0 ) < 0 and equality holds if and only if f(t) = 0 o.e. on T. 

Remark^ 1: The final condition establishes a geometric restriction on the vector A g 5R 2m . If we define the function 
2 : ^ ^ m by z(t) = $(t) -1 y(t) — $(t 0 )y 0 then A is restricted by the boundary condition \ T z(t f ) < 0. Moreover A is 

restricted by the whole trajectory, i.e. A T z(t) < 0 (( g T) and, upon differentiating, we obtain A T z(t) = — | q(t) \ f (i) a.e. on T. 
This shows that A is further restricted by the tangent vector to z(t) i.e. A T i(t) < 0 a.e. on T. Finally we observe that A T z{t) is 
monotone decreasing on T. All of these conditions place restrictions on the shape and location of an optimal trajectory. 

Remark 2: If the real number t 0 whose existence is asserted by this theorem is zero, then i) cannot hold; therefore (25, 
26) of ii) must hold where (26) requires that f(t) = 1 a.e. on T. A solution where t a = 0 is called an abnormal solution. 

On the other hand, a solution in which t 0 > 0 is called a normal solution. For normal solutions the constant l o /0 can be 
absorbed by the function q [22], Effectively, this is equivalent to setting l a = 0 in the proceeding theorem. For normal solutions 
the equation | q(t) |= 1 either has at most finitely many solutions on T and (26) is replaced by 


m x_f 0 . 1 ?(<)!< 1 
m " 1 1 . 1 *(*) l> 1 

(27) 

a.e. on T, or else 


1 9(01=1 (t€T) 

(28) 

and f is any real valued measurable function satisfying 


0 < /(<) < 1 

(29) 

a.e. on T such that (24, 25, 20) are satisfied. 
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3.3 Singular Solutions 

We shall say that a normal solution of the problem defined by (18 - 21) is singular if (28) holds, otherwise we say that it is 
nonsingular. Since (28) cannot hold for abnormal solutions, we shall also call them nonsingular. It follows from the preceeding 
theorem that (27) must hold for normal nonsingular solutions or else u(t) = 0 a.e. on T. For this reason and the continuity of 
q we can say that a normal nonsingular solution consists of thrusting intervals and coasting intervals separated by points on T 
called switches. It is sometimes convenient to refer to the switching function as the real valued function s on T defined by 

s(t) =| q(t) | -1. (30) 

The thrusting intervals are determined by the condition that the switching function is positive; similarly a negative value 
determines a coasting interval. A switch is a zero of the function s such that every open interval containing it also contains 
points where s is positive and points where s is negative. The preceeding theorem asserts that there are at most finitely many 
switches. The number of switches for a given interval T is not known, even for the linearized problem of a spacecraft near a 
circular orbit [28]. Since Theorem 3.1 shows that it is impossible to have s(t) = 0 on a subset of T of positive measure and also 
s(t) ^ 0 on another subset of T of positive measure, we do not need to define singular solutions on a subset of T. This property 
is lost, however, if the differential equation (18) is replaced by certain nonlinear differential equations such as those representing 
spacecraft trajectories in a Newtonian gravitational field, as discovered by Robbins [30] who found singular and nonsingular 
regimes on the same interval. If a solution is singular then (29) is valid a.e. on T, but the form of the function / is not specified. 
For this reason one might suspect that singular solutions to a boundary value problem are not unique. We shall show that this 
is the case for a certain large class of singular solutions. A singular solution is called an intermediate thrust solution if there 
exists a measurable subset S of T of positive measure such that 

0 < f(t) < 1 ( 31 ) 

for each t e S. We shall show that the intermediate thrust solutions are degenerate in the sense of the following theorem. This 
type of degeneracy was first discovered by La Salle [31] in the problem of time optimal control of linear systems. The theorem 
generalizes previous results based on linearized equations about a satellite in circular orbit [32]. An earlier observation of the 
degeneracy of singular solutions near a circular orbit was rrfctde by Marec [19]. t 

Theorem 3.2 Suppose that the assumptions of Theorem 3.1 are satisfied. If the problem of minimizing J[u] = f t * \ u(t) | dt 
subject to (18 - 20) has an intermediate thrust solution, then it has infinitely many intermediate thrust solutions. 

The proof can be found in [22]. 

3.4 The Two Point Boundary Value Problem 

In order to apply Theorem 3.1 one must solve a two-point boundary-value problem. A value of the element A 6 must be 
found such that the trajectory (24) satisfies the terminal condition (20). Specifically A defines the primer vector function q , 
which for normal nonsingular solutions defines the control function u through (25,27), establishing the trajectory through (24). 
We see therefore that the terminus y(tf) of the trajectory is a function of A G The terminal condition (20) can therefore 

be viewed as the solution of 2m nonlinear equations in the 2m components of A. 

Newton’s method is a well known computational method for the solution of several nonlinear equations in several variables, 
but reliable numerical convergence to a root requires knowledge of the approximate location of the root. We use a method 
of approximating the location of the root through the solution of the related impulsive problem [21]. The related impulsive 
problem provides a very good approximation if the number (3 in (18) is sufficiently large. The mathematical justification that 
the unbounded thrust problem can be viewed as the limit of the bounded thrust problem as ft tends to infinity is due to Neustadt 
[24]. The remarkable accuracy of the root for large values of ft is also partly due to the fact that the form of the primer vector 
function is found to be identical for both problems! Compare q in Theorem 3.1 with its corresponding definition in section 
2.2.1 Having solved the boundary-value problem for a large value of one can then successively lower the value of /?, solving 
the boundary value problem at each step, until the solution is reached for the designated value of /?. The whole sequence 
can usually be performed with only a microcomputer. The effectiveness of approximating a bounded thrust problem by an 
impulsive one and then lowering the bound in steps to solve the boundary value problem was demonstrated by Handelsman in 
1966 [33]. Starting iteratives for solution of the linear bounded thrust minimization problem can be obtained from solution of 
the related linear impulsive problem. A very attractive short cut is available, however, for many of the most useful problems 
in rendezvous and orbital maneuvers of spacecraft. For many problems in which </ - | v 0 |, and | v/ | are not too large, 

nonsingular bounded thrust solutions consist of at most two short duration thrusting intervals separated by a relatively large 
coasting interval. Starting iteratives are obtained for these problems from the related impulsive problem having at most two 
impulses at the end points, an initial increment Av 0 € & m at U and a terminal increment Av f 6 at t f . Computation of 
A v 0 and A Vf requires only the solution of the following set of 2m linear equations as determined by (14) and in [21] 

R(t 0 )Av 0 + R(tf)Av f = Qitfr'vf - *(<o)-V ( 32 ) 
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3.4.1 Algorithm for Computation of Optimal Trajectories Having Two Short Thrusting Intervals 

We now outline an algorithm for computation of certain two-point boundary-value problems associated with the solution of the 
optimal linear bounded thrust problem. The method is applicable to problems in which y 0 ,yj, and T satisfying the following: 

1. A nonsingular solution over the interval T exists connecting the points y a and y, and these points are not sufficiently 
close to points in which the only solutions are singular. It has been shown that, in some situations, the only solutions 
connecting Vo and y s are singular [32], and in these situations the algorithm is not applicable. Since the definition of q in 
Section 2.2.1 and (28) show that singular solutions are characterized by 

\ T R(t)R(t) T \ = 1 (t g T) (33) 

it is frequently possible to determine which elements A € Si 2 "* determine singular solutions and avoid those cases. 

2. The magnitudes | t»o |, | vj |, and tf — t 0 are small enough that nonsingular solutions consist of relatively short initial and 
terminal thrusting intervals separated by a relatively long coasting interval. It has been shown for the linear problem of 
optimal maneuvers near circular orbit that given any nontrivial interval T there are boundary values where more than 
two thrusting intervals occur; also more can occur for large values of t,-t 0 [8], It is conjectured for that problem that 
no more than four thrusting intervals can occur during one period of the circular orbit [28], For problems such as these 
the algorithm is not applicable, and the boundary-value problem is more difficult. 

The algorithm is outlined below. The process should generally begin using a much larger value of the number /? than is 
given in the problem. 

Step 1. The linear system (32) is solved for Av 0 and A v/. In some cases the 2m x 2m matrix (72(< 0 ), R(tf)) is singular. In 
these cases the value of tf should be changed by a small amount so that the resulting matrix is nonsingular. 

Step 2. The length of the initial and terminal thrusting intervals A t 0 and A tf are approximated by 

A to =| A v 0 | //?, A tj =| A Vf | //?. 

The switches t ao and b/ that define the coasting interval are 


— to + A t oy t a j = tf — A tj. 

Step 3. The values t Q and t f are replaced in (32) respectively by t ao and t a f t and steps 1 and 2 are repeated. If the 
magnitude of the difference in sucessive values of t ao and t s f satisfies a specified tolerance we proceed to step 4, otherwise steps 
1 and 2 are repeated. If the tolerance is not satisfied after a specified number of loops, the value of ft can be increased and the 
process can begin again at step 1. 

Step 4. The current values of Av c and A Vf are used to calculate the primer vector at the current values of t ao and t 3 f 
respectively. 

q(t*o) = -Av 0 / | Av 0 |, q(t a f) = - Av f / | Av f | . 

These equations hold for the related impulsive problem defined on the interval \t a0i t»j] [21]. 

Step 5. Utilizing the definition of q in Section 2.2.1 we can find the element A by solving the linear equations 


(R{t a o),R{t a ,)) T A 



If the matrix of coefficients is singular a very small adjustment in t a f is made t and the system is then solved for A. 

Step 6. Knowing A, the primer vector and its derivative are determined. The control function u can be calculated from 
(25) where /(<) is zero on the interval [t ao ,t a f] and /(<) = 1 otherwise. The vector y(tf)j& obtained through (11). If numerical 
integration of the differential equations describing y and q is preferred, then the initial values of the primer and its derivative 
are obtained from 


q(to) = R(tof A, q'(t 0 ) = R\t 0 ) T A. 

Step 7. The above yields y(tf). If this is sufficiently close to a root of the equation y(tf) ~y f = 0, then Newton’s method 
with iteration on A can be applied to determine a root. Standard commercial packages are available using shooting methods for 
numerical determination of the root, as well as performing numerical integration of the differential equations for q and y if this 
approach is preferred, and solving the systems of linear equations used herein [34]. After determination of the root, the value of 
/? is lowered, and a root is again found using the preceeding solution for starting iteratives. This process is repeated until /? is 
lowered to the correct value. 
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3.4.2 Spacecraft Maneuvers Near Circular Orbit 

We present some of the results for a 90 minute circular satellite orbit (w = 2*/3600 rad./sec.), a spacecraft mass of 3400kg, and 
a flight time of 600 sec. (f 0 = 0,t/ = 600). Although the maximum thrust of the spacecraft is 267 A, we assumed higher values 
initially, solved the boundary value problem, and successively lowered the thrust to its actual value. 

A typical example of the use of this algorithm is presented in figures 1 to 8. The transition from near impulsive to a low 
bounded thrust is shown graphically on the switching function curves where the portion above zero corresponds to thrusting 
and the portion of the curve below zero corresponds to coasting in figures 1 to 4. The trajectories required for fuel optimal 
rendezvous corresponding to these different thrust values are presented in figures 5 to 8. This particular calculation can be 
started at either 4000 A thrust where the maneuver requires thrusting for only 3% of the time which closely approximates the 
impulsive solution, or at 1000 A where the maneuver requires thrusting for 11% of the time. In either case the estimates of the 
starting values of Lawden’s primer vectors are close enough to obtain a converged solution. The thrust is then systematically 
lowered after each converged solution until a thrust of 267A is reached that requires thrusting for 65% of the maneuvering time. 

In general the starting primer vectors change very slowly as the maximum thrust is changed as long as the magnitude of 
thrust is large enough for the thrusting time to be less than 10% of the flight time. When the impulsive calculation predicts 
thrusting times greater than 20% of the flight time, then the changes in the magnitude of thrust between converged solutions 
produce rapidly increasing changes in the starting primer vectors. As the thrusting time approaches 50% of the flight time, 
progress becomes difficult unless the thrust magnitude is limited to small changes between converged solutions 

The calculations were made on an IBM PS/2 Model 50 computer which is on a Token Ring Network with IBM Personal 
Computer DOS Version 3.30. The programs are written and compiled with Microsoft(R) Fortran Optimizing Compiler Version 
4.10. The calculations are in single precision. 


4 A Fundamental Matrix Solution for Spacecraft Maneuvers Near a 
General Keplerian Orbit 

4.1 Lawden’s Integral and Some Previous Work for Noncircular Keplerian Orbits 

As early as 1954 Lawden [35] introduced a change of the variable from time to the true anamaly in describing the equations of 
motion of a rocket under the action of an inverse square law of force, and in order to evaluate his primer vector during a nu 
thrust (Keplerian) interval, he introduced the integral 


m 


-f 

Je 0 


dj 

(1 + e cos £) 2 sin 2 $ 


(34) 


where 0 denotes the true anamaly at time t and 0 O denotes its initial value at time t 0 . In 1963 this integral appears again in his 
book [25] in connection with the solution of his equations (5.30) - (5.32) which determine the transformed primer vector along 
a Keplerian interval. 

These same equations, with only minor modifications, appeared at approximately the same time in the work of De Vries 
[9] to describe the relative motion of two nearby points in elliptical orbits. Except for changes in the coordinate systems, a 
nonhomogeneous version of these equations was used by Tschauner and Hempel [10] to describe the rendezvous of a spacecraft 
with a target in elliptical orbit. These equations were used again by Tschauner [14] and were investigated through a c ange o 
variable from the true anamaly to the eccentric anamaly. Restricted or generalized forms of these equations were employed by 
Shulman and Scott [11], Euler and Shulman [13], and Euler [12] for rendezvous of a spacecraft with an object in elliptical orbit. 
The approach of Tschauner and Hempel, as used by Weiss [15] in 1981 was found to be effective in constructing two-impulse 
solutions to rendezvous problems involving objects in elliptical orbits of high eccentricity [16]. In all these stu ies, t e so ution 
of Lawden’s equations was not investigated through the use of the integral 1(0) of Eq. (34). 

Eckel’s paper [36] of 1982 returns to the integral 1(0) with Lawden’s equations to determine the primer vector and solve 
the problem of optimal impulsive transfer between noncoplanar elliptical orbits. Recently Carter and Humi [17], and Carter [29] 
used this integral in solving Lawden’s equations to determine both the primer vector and the structure of an optimal ren ezvous 
of a spacecraft with an object near a point in general Keplerian orbit. ? 

The previously mentioned studies demonstrate a variety of applications and interpretations of Lawden s equations. They 
also indicate considerable variation in the form of the solution of these equations. 

We prefer those forms which include the integral 1(0) because it appears naturally in the most straightforward approaches 
to solving Lawden’s differential equation. This integral is singular, however, at values where the true anamaly is a mu tip e 
of t. Even though these are removable singularities, they may lead to computational instabilities in the solution. Lspeci y 
bothersome is the fact that computational problems can occur where the true anamaly is near zero. These pro lems are avoi e 
in the work of Tschauner and Hempel [10,14] and others [11-13, 15] who use a form of solution that does not involve 1(0), but 
their work is confined to elliptical orbits. 


550 



*.*0 


SwiUhtag „ 
p*miu. ,0 * 



•1X0 1 1 

0.00 M0X0 400X0 

Tim* ot flight (mc) 



ri|w« 1. IwlUhiag Puacllaa »». night Tim* tot Thruat ot 4000N, 
X t {0) a* 0,X,(0) - 1000m, Jf»(«00) - 0, X,(«00) - 0. 


P " U f#r Tfc '«» ■» 4000N, 
*»(°> - 0,X,(0) - 1000m, X,{«00) . 0, *,(«#) - 0. 




rtgur* J. Switching Puactloa vs. Plight Tim* tot Thrust ot 1000N, 

*»(*) - 0, X»(0) - 1000m, X, (BOO) - O.X, (100) - 0. 


Flfur* I. Spacecraft Plight Path tor Thrust ot 1000N 
X,(0) - 0, JV|(0) - 1000m, ATi(OOO) - 0,X,(600) as 0. 


SwIUkUg 

PuaclUa 




Figure 9. SwiUhlag Puadloa va. night Tima tot Thruat ot SOON. 
X|(0) - 0, *,(0) - 1000m, X|(000) - 0, X,(«O0) - 0. 


Figvra 7. Spacecraft Plight Path (hr Thruat ot WON, 
*i(0) - 0,Jf,(0) - 1000m, X,(S00) • 0,X,{«00) - 0, 




Pig or* 4. Switching PhacUmt **. PUgkt Tim* for Thruat ot SI7M, 
X,(0) - O.X,(0) a. 1000m, X, (600) - 0,X,{«00) « 0. 


Figur* I. Spacecraft Plight Path lor Thrust ot M7N, 
X|(0) - 0. X,(0) at 1000m, X, (600) as 0,X,(600) - 0. 


551 


ORIGINAL PAGE IS 
OF POOR QUALITY 








The approach presented in this paper is to modify the original form by replacing 1(9) by a related integral J(9), thereby 
removing all singularities and computational instabilities. The resulting solution, in terms of J(9), is identical for hyperbolic, 
parabolic, or noncircular elliptic orbits, but the particular case determines the nature of the closed-form evaluation of J(6). 

4.2 Transformed Equations of a Spacecraft near Keplerian Orbit 

We consider a rotating coordinate frame centered at a point moving in a Keplerian orbit about a central attractive body. The 
positive x 2 axis is directed away from the central body, the positive x, axis is perpendicular to it and opposes the direction 
of the motion, and the x 3 axis completes a right handed coordinate system. We consider the equations of motion of a point 
mass spacecraft relative to this coordinate system in which the Newtonian gravitational force function has been linearized about 
the point in Keplerian orbit. The independent variable is the true anamaly 0 defined on the closed interval 8 0 < 8 < 9, 
which we denote by 6. All vectors are assumed to be elements of three dimensional Euclidean space The position vector 
x(9) = (xi(0),x 2 (0),X3(0)) of the spacecraft in this coordinate system is transformed to the vector z(9) = (zi(9),z 2 (9),z 3 (9)) 


by the equation 


z(0) = r{9)x(0) 


(35) 


where 

r(0) = 1 + ecos(9) 


(36) 


and e denotes the eccentricity of the Keplerian orbit. 

This development which is presented in detail in previous work [17] results in the following transformed equations of the 
spacecraft 


*"(*) = 2*7(0) + *iW) 


*2(») = - 2*5 («) + *2(9) 

z"(9) = -z 3 (9) + q 3 (9) 

where the prime indicates differentiation with respect to 8, and the vector a(9) = (<»i(0),a 2 (0),a 3 (0)) is given by 

a(9) = m*(9) 


(37) 

(38) 


where , 

m = k/r(0) 3 . 


(39) 


In this expression the positive constant k is I 6 T m //i 4 m where L is the magnitude of the constant angular momentum of 
the object in Keplerian orbit divided by its mass, T m is the maximum magnitude of the thrust of the spacecraft, p is the product 
of the universal gravitational constant and the mass of the central body of attraction, and m is the mass of the spacecraft. The 
vector u(0) = (t*i(0), u 2 (0), u 3 (0)) represents the normalized thrust of the spacecraft and is subject to the constraint 


I *(•) |< 1. 


(40) 


Since this investigation is restricted to linear equations, we shall assume that the mass m is constant. Previous investigations, 
however, have taken into account the change in mass of the spacecraft as propellant is consumed [29,18]. 

Equations (37) are essentially the equations of Tschauner and Hempel [10], and their homogeneous form represents essen- 
tially the equations of De Vries [9] and Lawden [25]. 

Here the class of admissible control functions is the set of all Lebesque measurable vector valued functions that satisfy (40) 
a.e. on ©. The optimal rendezvous problem associated with a point in Keplerian orbit is defined as the determination of an 
admissible control function t* that minimizes the cost function 



subject to the conditions (37-40) which are valid a.e. on 0 and the end conditions 


z(9 0 ) = zo, z'(9 0 ) = vo, 

z(9 s ) = zj y z'(9f) = v f ( 42 ) 

where the vectors zo and z / define the initial and terminal values of the transformed position z, and vo and v / define the initial 

and terminal values of its derivative. . 

Theorem 3.1 is not directly applicable to this problem because (38) and (4.2) do not conform to (18) and (21) respective y. 
The theorem can be generalized to include this problem, however, and the results are very similar if the primer vector q(0) is 
replaced by the transformed primer Q(9) = q(9)/r(9) as in previous work. [17,29] 
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If we define the state vector function y : © 
similar to (18) where 0(8) replaces 0 and 


** b y y{8) = (z(0) T , zl(0) T ) T , then equations (37-39) can be put in a form 


m = 


0 0 0 1 0 0 * 

0 0 0 0 1 0 

0 0 0 0 0 1 

0 0 0 0 2 0 

0 3/r(0) 0-200 

0 0 -1 0 0 0 


(43) 


Here 8 replaces t and © replaces T. By defining j, 0 = (zf,v 0 T ) T and y, = (zJ,vj) T the boundary conditions (41) satisfy (19, 

In order to obtain a fundamental matrix solution *(8) associated with this problem, we solve the homogeneous form of 
the system (37), that is, we find the complete solution where a(9) is identically zero. This solution was found by Lawden [35] in 
terms of the integral 1(B). Using a similar form of solution [17], we obtain 


Zl (0) = -61 r(9) 2 - b 2 [r(8) 2 I(8) + cotS ] - 6 3 stn(?[l + r(0)] + 64 


z 2 {8) = er(fl)stn^[6 1 + b 2 I(0)] - b 3 r(8)cos0 


z 3 (B) = bscosB + besinB (44) 

where 61 , b 2i b 3 , b 4 , bs,b 6 are arbitrary constants of integration. 

These equations can be used to define the relative motion of two nearby points in non circular Keplerian orbits by trans- 
forming from z(B) to x(B) using (35), and is a generalization of the work of De Vries. [9]. 


4.3 New Form of the Rendezvous Equations and the Resulting Fundamental Matrix 
Solution 

The Equation (43) has removable singularities at integer multiples of x which we denote by #r. These singularities appear in 
the expressions 1(B) and cotB. Computation is troublesome at or near these singularities. These computational instabilities can 
be removed. If we integrate (34) by parts, it becomes 


m = 2eJ(8)-^ + c ( 45 ) 

which holds except at the singularities of 1(0) where c is a constant of integration and 



(46) 


It is observed that the integral J(0) has no singularities. It follows from the continuity of J(0) that 

2eJ(„x) + C = tf Hm [/(*) + (47) 

so that the singularities in (43) are removed. It follows that there are no singularities if solutions of (37) are stated in terms of 
J(B) instead of 1(B). [18] With these changes the solution (43) of the unpowered equations of motion becomes 


zi(B) = -r(B) 2 [b! + 2heJ(9)] - ft 3 [l -f r(B)]sinB + 


Z2(B) — r(^)s*n^[ 6 1 e + 2b2t 2 J(B)] — cosB[~^ -|- 63 r(^)] 

*$(0) = b^cosB -h besinB. ( 4 g) 

These equations are valid for all noncircular Keplerian orbits. The integral J(8) can be evaluated in closed form, and the 
particular form is determined from the type of orbit. Equations (47) are more useful than equations (43). If we differentiate 
(47) we obtain the vector-valued function z' : © S? 3 . The pair (z(8) T ,zi(9) T ) T defines the state vector y(9). The equation 

z(0) from (47) and its derivative z'(0) together are represented by 


y(0) = $(0)b (0e&) 


(49) 
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where b = (bub 2 ,b 3 ,b t ,h,be) T is a constant vector in » 6 and $(0) is a fundamental matrix solution associated with A(6). It 
follows from (47) that 


*(*) = 


-r(«) 3 
er(# ) sin 9 
0 

3er(0) ain 9 
er(«) coa 9 — e 3 ain 
0 


3©r(«)ain# J(9) - co*9/r(9 ) 

0 

4er(*)ain*J(«) - 3 coa 9fr(9) 

^ »in#(l+2eco»<) 

3e(r(0) coa 9 — e ain 3 0) J(0) + r(9 ) 3 


“(1 + »•(*)) 

— r(0) coa 9 
0 

e — 3r(0) coa 9 
(r(0) + c coa 9 ) mind 
0 


10 0 
0 0 0 

0 coa 9 ain# 

0 0 0 

0 0 0 

0 — ain 9 coa 9 


(50) 


Our goal has been to determine this matrix function. If the problem is restricted to a plane, we delete the third and sixth row and 
the fifth and sixth column. The resulting 4x 4 matrix function is a fundamental matrix solution for the planer problem. There 
are several approaches 3 to the problem of inverting a fundamental matrix solution and applying a generalization of the preceeding 
theory. These should result in new methods of computing either impulsive or bounded thrust trajectories of a spacecraft near 
Keplerian orbit. Further results in this area are forthcoming. 


4,4 Closed-Form Evaluation of The New Integral 

We show here that the integral J(9) can be evaluated easily if we transform from the true anamaly 8 to the eccentric anamaly 
E for elliptical orbits or its analog H for hyperbolic orbits. 


4.4.1 Elliptical Orbits 

For orbits in which 0 < e < 1 we have the relationship between the eccentric anomaly E and true anomaly 0 given by 

. COS E C / c i \ 

COS 0 = “ (51) 

1 — e cos E 

where sin 6 and sin E always have the same algebraic sign. Changing the variable to E in (45) establishes the much simpler 
integral £ 

7(0) = ( l-e 2 ) _s/2 / (1 - ecos£)(cos£ - e)d£ (52) 

Je 0 

where E 0 is the eccentric anomaly at 0 O . This integral is easily evaluated using elementary methods to obtain 


J(0) = -(1 - e 2 ) _5/2 [yU - (1 + e 2 )sin E + | sin E cos E + cj 
where C is an arbitrary constant. 


(53) 


4.4.2 Hyperbolic Orbits 

In a similar way we can evaluate J(0) for orbits in which e > 1. We introduce the analog of the eccentric anomaly H by the 
relationship 

cos 0 = c ~ coshg . (54) 

e cosh H — 1 

where sin 0 and sinh H always have the same algebraic sign. With this substitution the integral defined by (45) becomes 


r H 

J(9) = (e 2 — i)~ 5 / 2 / (e cosh £ — l)(e — cosh £)d£ 

Jh 0 

where H 0 denotes the value of H at 0 O . This integral also is easily evaluated. The result is 

J(0) = -(e 2 - l)- s/2 - (1 + e 2 ) sinh H + cosh H + c] 

where again C denotes an arbitrary constant. 


(55) 


(56) 


4.4.3 Parabolic Orbits 

For the case in which e = 1, the integral J(&) can be evaluated directly using the identity cos 6 = 2 cos 2 f — 1. The result is 

J(0) = \ + C ( 57 ) 

where C is again an arbitrary constant. Since this expression is only defined on the region r(0) > O(i.e.cos0 > —1), it has no 
singularities. 

3 This matrix was inverted by Prof. Mayer Humi of the Mathematial Sciences Department at Worchester Polytechnic Institute using 
MAXYSMA. The author has also inverted this matrix using the adjoint system to that defined by (42). Numerical inversion is also feasible. 
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